1 Set Up

#Custom function that finds factors with only 1 level
#Source: https://stackoverflow.com/questions/44200195/how-to-debug-contrasts-can-be-applied-only-to-factors-with-2-or-more-levels-er

debug_contr_error <- function (dat, subset_vec = NULL) {
  if (!is.null(subset_vec)) {
    ## step 0
    if (mode(subset_vec) == "logical") {
      if (length(subset_vec) != nrow(dat)) {
        stop("'logical' `subset_vec` provided but length does not match `nrow(dat)`")
        }
      subset_log_vec <- subset_vec
      } else if (mode(subset_vec) == "numeric") {
      ## check range
      ran <- range(subset_vec)
      if (ran[1] < 1 || ran[2] > nrow(dat)) {
        stop("'numeric' `subset_vec` provided but values are out of bound")
        } else {
        subset_log_vec <- logical(nrow(dat))
        subset_log_vec[as.integer(subset_vec)] <- TRUE
        } 
      } else {
      stop("`subset_vec` must be either 'logical' or 'numeric'")
      }
    dat <- base::subset(dat, subset = subset_log_vec)
    } else {
    ## step 1
    dat <- stats::na.omit(dat)
    }
  if (nrow(dat) == 0L) warning("no complete cases")
  ## step 2
  var_mode <- sapply(dat, mode)
  if (any(var_mode %in% c("complex", "raw"))) stop("complex or raw not allowed!")
  var_class <- sapply(dat, class)
  if (any(var_mode[var_class == "AsIs"] %in% c("logical", "character"))) {
    stop("matrix variables with 'AsIs' class must be 'numeric'")
    }
  ind1 <- which(var_mode %in% c("logical", "character"))
  dat[ind1] <- lapply(dat[ind1], as.factor)
  ## step 3
  fctr <- which(sapply(dat, is.factor))
  if (length(fctr) == 0L) warning("no factor variables to summary")
  ind2 <- if (length(ind1) > 0L) fctr[-ind1] else fctr
  dat[ind2] <- lapply(dat[ind2], base::droplevels.factor)
  ## step 4
  lev <- lapply(dat[fctr], base::levels.default)
  nl <- lengths(lev)
  ## return
  list(nlevels = nl, levels = lev)
  }

#Variables to put into model

  • Dose (counts)
  • Dose (mass)
  • Dose (volume)
  • Size (continuous)
  • Size (binned)
  • Polymer type
  • Shape
  • Organism Group (this is being used as a temporary substitute for body size)
  • Life Stage
  • Level of Biological Organization
  • Exposure Duration
  • Acute/Chronic (this values are only complete for fish, molluscs, crustaceans and algae)
  • Charge (negative/positive - categorical)
  • Zeta potential
  • Particle Source (commercial, generated in-house, or mmodified commercial particles e.g., milling)

Response Variables * Effect (Y/N) - Binary Data * Effect (Y/N) - Binary Data x Effect Score (binned organism level effects)

#Continous Variable Distributions

There are several continuous variables that we want to feed into our model. Before doing so, we’re going to check the distribution to see if any of them are skewed and need to be transformed.

Due to skewed data, the following categories are log10 transformed before modeling:

  • Dose (counts)
  • Dose (mass)
  • Dose (volume)
  • Size (continuous)
  • Exposure Duration

2 Model: Kitchen Sink

All Independent Variables, Response Variable: Effect (Y/N)

Full Model

## $nlevels
##         shape_f           org_f          life_f           bio_f acute.chronic_f 
##               3               3               3               5               2 
## 
## $levels
## $levels$shape_f
## [1] "Fiber"    "Fragment" "Sphere"  
## 
## $levels$org_f
## [1] "Crustacea" "Fish"      "Mollusca" 
## 
## $levels$life_f
## [1] "Adult"    "Early"    "Juvenile"
## 
## $levels$bio_f
## [1] "Cell"       "Organism"   "Population" "Subcell"    "Tissue"    
## 
## $levels$acute.chronic_f
## [1] "Acute"   "Chronic"
## 
## Call:
## lm(formula = effect_10 ~ ., data = aoc_setup_select_1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6552 -0.3369 -0.1877  0.4716  1.1018 
## 
## Coefficients:
##                                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)                              1.051687   0.150037   7.010 2.87e-12
## log.dose.particles.mL.master            -0.112899   0.125461  -0.900  0.36825
## log.dose.mg.L.master                     0.051598   0.008582   6.012 2.03e-09
## log.dose.um3.mL.master                   0.093169   0.125761   0.741  0.45884
## log.size.length.um.used.for.conversions -0.361561   0.376854  -0.959  0.33742
## shape_fFragment                         -0.261101   0.046653  -5.597 2.36e-08
## shape_fSphere                           -0.271240   0.101701  -2.667  0.00769
## org_fFish                               -0.071743   0.023682  -3.029  0.00247
## org_fMollusca                           -0.165614   0.026550  -6.238 4.98e-10
## life_fEarly                             -0.054630   0.022028  -2.480  0.01319
## life_fJuvenile                          -0.110397   0.023656  -4.667 3.18e-06
## bio_fOrganism                           -0.341509   0.047869  -7.134 1.18e-12
## bio_fPopulation                         -0.616304   0.095914  -6.426 1.50e-10
## bio_fSubcell                            -0.124617   0.045353  -2.748  0.00603
## bio_fTissue                             -0.167895   0.055423  -3.029  0.00247
## log.exposure.duration.d                  0.058345   0.014896   3.917 9.15e-05
## acute.chronic_fChronic                   0.051101   0.023834   2.144  0.03210
##                                            
## (Intercept)                             ***
## log.dose.particles.mL.master               
## log.dose.mg.L.master                    ***
## log.dose.um3.mL.master                     
## log.size.length.um.used.for.conversions    
## shape_fFragment                         ***
## shape_fSphere                           ** 
## org_fFish                               ** 
## org_fMollusca                           ***
## life_fEarly                             *  
## life_fJuvenile                          ***
## bio_fOrganism                           ***
## bio_fPopulation                         ***
## bio_fSubcell                            ** 
## bio_fTissue                             ** 
## log.exposure.duration.d                 ***
## acute.chronic_fChronic                  *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4382 on 3385 degrees of freedom
## Multiple R-squared:  0.1104, Adjusted R-squared:  0.1062 
## F-statistic: 26.25 on 16 and 3385 DF,  p-value: < 2.2e-16

Stepwise Model - Both Directions

## 
## Call:
## lm(formula = effect_10 ~ log.dose.particles.mL.master + log.dose.mg.L.master + 
##     log.size.length.um.used.for.conversions + shape_f + org_f + 
##     life_f + bio_f + log.exposure.duration.d + acute.chronic_f, 
##     data = aoc_setup_select_1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6616 -0.3384 -0.1883  0.4722  1.1004 
## 
## Coefficients:
##                                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)                              0.961328   0.087373  11.003  < 2e-16
## log.dose.particles.mL.master            -0.020160   0.008385  -2.404  0.01626
## log.dose.mg.L.master                     0.051909   0.008571   6.056 1.55e-09
## log.size.length.um.used.for.conversions -0.083029   0.025820  -3.216  0.00131
## shape_fFragment                         -0.268834   0.045467  -5.913 3.70e-09
## shape_fSphere                           -0.205081   0.048657  -4.215 2.57e-05
## org_fFish                               -0.071653   0.023680  -3.026  0.00250
## org_fMollusca                           -0.164836   0.026527  -6.214 5.80e-10
## life_fEarly                             -0.053066   0.021925  -2.420  0.01556
## life_fJuvenile                          -0.109698   0.023636  -4.641 3.59e-06
## bio_fOrganism                           -0.341726   0.047865  -7.139 1.14e-12
## bio_fPopulation                         -0.616644   0.095906  -6.430 1.46e-10
## bio_fSubcell                            -0.124734   0.045350  -2.750  0.00598
## bio_fTissue                             -0.168007   0.055420  -3.032  0.00245
## log.exposure.duration.d                  0.058166   0.014893   3.905 9.59e-05
## acute.chronic_fChronic                   0.050621   0.023824   2.125  0.03368
##                                            
## (Intercept)                             ***
## log.dose.particles.mL.master            *  
## log.dose.mg.L.master                    ***
## log.size.length.um.used.for.conversions ** 
## shape_fFragment                         ***
## shape_fSphere                           ***
## org_fFish                               ** 
## org_fMollusca                           ***
## life_fEarly                             *  
## life_fJuvenile                          ***
## bio_fOrganism                           ***
## bio_fPopulation                         ***
## bio_fSubcell                            ** 
## bio_fTissue                             ** 
## log.exposure.duration.d                 ***
## acute.chronic_fChronic                  *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4381 on 3386 degrees of freedom
## Multiple R-squared:  0.1102, Adjusted R-squared:  0.1063 
## F-statistic: 27.96 on 15 and 3386 DF,  p-value: < 2.2e-16

3 Model: Crustacean Fitness, Organism Level Endpoints Only

3.1 Dose (mass)

Full Model

## 
## Call:
## glm(formula = effect_10 ~ logdose.mg.L.master * size.length.um.used.for.conversions, 
##     family = "binomial", data = m1_crust, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8447  -0.7063  -0.6178  -0.5142   2.0485  
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                             -1.425e+00  7.388e-02
## logdose.mg.L.master                                      1.807e-01  4.837e-02
## size.length.um.used.for.conversions                     -3.807e-04  1.785e-04
## logdose.mg.L.master:size.length.um.used.for.conversions  1.057e-04  4.555e-05
##                                                         z value Pr(>|z|)    
## (Intercept)                                             -19.288  < 2e-16 ***
## logdose.mg.L.master                                       3.737 0.000187 ***
## size.length.um.used.for.conversions                      -2.133 0.032944 *  
## logdose.mg.L.master:size.length.um.used.for.conversions   2.320 0.020367 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1311.5  on 1310  degrees of freedom
## Residual deviance: 1272.9  on 1307  degrees of freedom
## AIC: 1280.9
## 
## Number of Fisher Scoring iterations: 4
#Plots

ggPredict(m1_crust_model_dose,colorn=100,jitter=FALSE, interactive = TRUE)

3.2 Dose (volume)

Full Model

## 
## Call:
## glm(formula = effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions, 
##     family = "binomial", data = m1_crust, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8116  -0.6866  -0.6282  -0.5292   2.0869  
## 
## Coefficients:
##                                                             Estimate Std. Error
## (Intercept)                                               -2.176e+00  2.574e-01
## logdose.um3.mL.master                                      1.262e-01  4.049e-02
## size.length.um.used.for.conversions                       -1.176e-03  4.464e-04
## logdose.um3.mL.master:size.length.um.used.for.conversions  1.259e-04  4.533e-05
##                                                           z value Pr(>|z|)    
## (Intercept)                                                -8.451  < 2e-16 ***
## logdose.um3.mL.master                                       3.118  0.00182 ** 
## size.length.um.used.for.conversions                        -2.635  0.00842 ** 
## logdose.um3.mL.master:size.length.um.used.for.conversions   2.777  0.00549 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1311.5  on 1310  degrees of freedom
## Residual deviance: 1277.6  on 1307  degrees of freedom
## AIC: 1285.6
## 
## Number of Fisher Scoring iterations: 4
#Plots
ggPredict(m1_crust_model_dose,
          colorn=200,
 #         jitter=TRUE, 
point = FALSE,
          interactive = TRUE, 
show.summary = TRUE,
#digits = 3, 
se = FALSE)
## 
## Call:
## glm(formula = effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions, 
##     family = "binomial", data = m1_crust, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8116  -0.6866  -0.6282  -0.5292   2.0869  
## 
## Coefficients:
##                                                             Estimate Std. Error
## (Intercept)                                               -2.176e+00  2.574e-01
## logdose.um3.mL.master                                      1.262e-01  4.049e-02
## size.length.um.used.for.conversions                       -1.176e-03  4.464e-04
## logdose.um3.mL.master:size.length.um.used.for.conversions  1.259e-04  4.533e-05
##                                                           z value Pr(>|z|)    
## (Intercept)                                                -8.451  < 2e-16 ***
## logdose.um3.mL.master                                       3.118  0.00182 ** 
## size.length.um.used.for.conversions                        -2.635  0.00842 ** 
## logdose.um3.mL.master:size.length.um.used.for.conversions   2.777  0.00549 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1311.5  on 1310  degrees of freedom
## Residual deviance: 1277.6  on 1307  degrees of freedom
## AIC: 1285.6
## 
## Number of Fisher Scoring iterations: 4
#particle couunt
#Select the columns that we want to feed into the model for simplicity
m1_crust_part <- aoc_setup %>% 
  filter(!effect_f == "NA") %>% 
  #filter(effect.metric == c("HONEC", "LOEC", "NOEC")) %>% 
  #filter(organism.group == "Crustacea") %>% 
  filter(!acute.chronic_f == "Unavailable") %>% 
  filter(bio.org == "organism") %>% 
  filter(!environment == "Terrestrial") %>% 
  filter(lvl1_f == "Fitness") %>% 
  filter(!size_f == "Not Reported") %>% 
  mutate(logdose.particles.mL.master = log10(dose.particles.mL.master)) %>% 
 # filter(acute.chronic_f == "Acute") %>% 
  dplyr::select(c(effect_10, logdose.particles.mL.master, size.length.um.used.for.conversions)) %>% 
  drop_na()

#Fit the full model
m1_crust_model_particles <- glm(effect_10 ~ logdose.particles.mL.master * size.length.um.used.for.conversions, data = m1_crust_part, na.action = "na.exclude", family = "binomial")

summary(m1_crust_model_particles)
## 
## Call:
## glm(formula = effect_10 ~ logdose.particles.mL.master * size.length.um.used.for.conversions, 
##     family = "binomial", data = m1_crust_part, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5942  -0.6724  -0.6214  -0.5551   2.0530  
## 
## Coefficients:
##                                                                   Estimate
## (Intercept)                                                     -1.808e+00
## logdose.particles.mL.master                                      7.178e-02
## size.length.um.used.for.conversions                              2.881e-04
## logdose.particles.mL.master:size.length.um.used.for.conversions  1.331e-04
##                                                                 Std. Error
## (Intercept)                                                      1.355e-01
## logdose.particles.mL.master                                      2.177e-02
## size.length.um.used.for.conversions                              6.364e-05
## logdose.particles.mL.master:size.length.um.used.for.conversions  4.499e-05
##                                                                 z value
## (Intercept)                                                     -13.342
## logdose.particles.mL.master                                       3.297
## size.length.um.used.for.conversions                               4.528
## logdose.particles.mL.master:size.length.um.used.for.conversions   2.959
##                                                                 Pr(>|z|)    
## (Intercept)                                                      < 2e-16 ***
## logdose.particles.mL.master                                     0.000978 ***
## size.length.um.used.for.conversions                             5.96e-06 ***
## logdose.particles.mL.master:size.length.um.used.for.conversions 0.003090 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1311.5  on 1310  degrees of freedom
## Residual deviance: 1280.0  on 1307  degrees of freedom
## AIC: 1288
## 
## Number of Fisher Scoring iterations: 4

plot

ggPredict(m1_crust_model_particles, colorn=1000, interactive = TRUE)

3.2.1 volume

#particle couunt
#Select the columns that we want to feed into the model for simplicity
m1_crust_volume <- aoc_setup %>% 
  filter(!effect_f == "NA") %>% 
  #filter(effect.metric == c("HONEC", "LOEC", "NOEC")) %>% 
  #filter(organism.group == "Crustacea") %>% 
  filter(!acute.chronic_f == "Unavailable") %>% 
  filter(bio.org == "organism") %>% 
  filter(!environment == "Terrestrial") %>% 
  filter(lvl1_f == "Fitness") %>% 
  filter(!size_f == "Not Reported") %>% 
  mutate(logdose.um3.mL.master = log10(dose.um3.mL.master)) %>% 
 # filter(acute.chronic_f == "Acute") %>% 
  dplyr::select(c(effect_10, logdose.um3.mL.master, size.length.um.used.for.conversions)) %>% 
  drop_na()

#Fit the full model
m1_crust_model_volume <- glm(effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions, data = m1_crust_volume, na.action = "na.exclude", family = "binomial")

summary(m1_crust_model_volume)
## 
## Call:
## glm(formula = effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions, 
##     family = "binomial", data = m1_crust_volume, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8116  -0.6866  -0.6282  -0.5292   2.0869  
## 
## Coefficients:
##                                                             Estimate Std. Error
## (Intercept)                                               -2.176e+00  2.574e-01
## logdose.um3.mL.master                                      1.262e-01  4.049e-02
## size.length.um.used.for.conversions                       -1.176e-03  4.464e-04
## logdose.um3.mL.master:size.length.um.used.for.conversions  1.259e-04  4.533e-05
##                                                           z value Pr(>|z|)    
## (Intercept)                                                -8.451  < 2e-16 ***
## logdose.um3.mL.master                                       3.118  0.00182 ** 
## size.length.um.used.for.conversions                        -2.635  0.00842 ** 
## logdose.um3.mL.master:size.length.um.used.for.conversions   2.777  0.00549 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1311.5  on 1310  degrees of freedom
## Residual deviance: 1277.6  on 1307  degrees of freedom
## AIC: 1285.6
## 
## Number of Fisher Scoring iterations: 4

plot

ggPredict(m1_crust_model_volume, colorn=500, interactive = TRUE)
#alternative method using ggplot

#get probabilities
m1 <- m1_crust %>% 
  mutate(prob = predict(m1_crust_model_dose,
                        type = "response"))
  
# plot
ggplot(m1, aes(x = logdose.mg.L.master, y = effect_10, color = size.length.um.used.for.conversions))+
  geom_jitter() +
  geom_smooth(method = 'glm',
              aes(x = logdose.mg.L.master, y = prob), #use predicted NOEC probability
              method.args = list(family = binomial(link = 'logit')),
              se = FALSE, color = 'red') #+

  # geom_smooth(method = 'glm',
  #             method.args = list(family = binomial(link = 'probit')),
  #             se = FALSE, color = 'green') +
  # stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE)
ggPredict(m1_crust_model_dose,colorn=100,jitter=FALSE, interactive = TRUE)

3.3 Dose (count)

Full Model

## 
## Call:
## glm(formula = effect_10 ~ logdose.particles.mL.master * size.length.um.used.for.conversions, 
##     family = "binomial", data = m1_crust, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5942  -0.6724  -0.6214  -0.5551   2.0530  
## 
## Coefficients:
##                                                                   Estimate
## (Intercept)                                                     -1.808e+00
## logdose.particles.mL.master                                      7.178e-02
## size.length.um.used.for.conversions                              2.881e-04
## logdose.particles.mL.master:size.length.um.used.for.conversions  1.331e-04
##                                                                 Std. Error
## (Intercept)                                                      1.355e-01
## logdose.particles.mL.master                                      2.177e-02
## size.length.um.used.for.conversions                              6.364e-05
## logdose.particles.mL.master:size.length.um.used.for.conversions  4.499e-05
##                                                                 z value
## (Intercept)                                                     -13.342
## logdose.particles.mL.master                                       3.297
## size.length.um.used.for.conversions                               4.528
## logdose.particles.mL.master:size.length.um.used.for.conversions   2.959
##                                                                 Pr(>|z|)    
## (Intercept)                                                      < 2e-16 ***
## logdose.particles.mL.master                                     0.000978 ***
## size.length.um.used.for.conversions                             5.96e-06 ***
## logdose.particles.mL.master:size.length.um.used.for.conversions 0.003090 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1311.5  on 1310  degrees of freedom
## Residual deviance: 1280.0  on 1307  degrees of freedom
## AIC: 1288
## 
## Number of Fisher Scoring iterations: 4
require(ggiraphExtra)
#Plots
ggPredict(m1_crust_model_dose,colorn=100,jitter=FALSE, interactive = TRUE)

3.4 Discrete predictor and dose

## 
## Call:
## glm(formula = effect_10 ~ logdose.mg.L.master * size_f, family = "binomial", 
##     data = m1_crust, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7085  -0.7062  -0.6167  -0.4255   2.3881  
## 
## Coefficients:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                           -1.57691    0.21094  -7.476 7.67e-14 ***
## logdose.mg.L.master                    0.40501    0.16344   2.478   0.0132 *  
## size_f100nm < 1µm                     -0.01864    0.37454  -0.050   0.9603    
## size_f1µm < 100µm                      0.17395    0.22858   0.761   0.4467    
## size_f100µm < 1mm                     -0.25876    0.35626  -0.726   0.4676    
## size_f1mm < 5mm                       -0.87511    0.64223  -1.363   0.1730    
## logdose.mg.L.master:size_f100nm < 1µm  0.45050    0.31359   1.437   0.1508    
## logdose.mg.L.master:size_f1µm < 100µm -0.24482    0.17319  -1.414   0.1575    
## logdose.mg.L.master:size_f100µm < 1mm -0.23277    0.22994  -1.012   0.3114    
## logdose.mg.L.master:size_f1mm < 5mm    0.04129    0.21904   0.189   0.8505    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1358.5  on 1347  degrees of freedom
## Residual deviance: 1303.8  on 1338  degrees of freedom
## AIC: 1323.8
## 
## Number of Fisher Scoring iterations: 5

Plotted below.

discrete <- ggpredict(m1_crust_model_discrete,terms = c("logdose.mg.L.master", "size_f"))
plot(discrete, facet = TRUE, add.data = TRUE, colors = "quadro", alpha = 0.2, dodge = 0.2)

     #title = "All Aquatic Organisms")

3.4.0.1 Acute only crustacea

## 
## Call:
## glm(formula = effect_10 ~ logdose.mg.L.master * size_f + logdose.particles.mL.master * 
##     size_f, family = "binomial", data = m1_crust_acute, na.action = "na.exclude")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4072  -0.8096  -0.6639  -0.2631   2.1764  
## 
## Coefficients:
##                                               Estimate Std. Error z value
## (Intercept)                                    -30.504      9.991  -3.053
## logdose.mg.L.master                             -2.914      1.075  -2.711
## size_f100nm < 1µm                               29.153     10.270   2.839
## size_f1µm < 100µm                               28.853      9.996   2.886
## size_f100µm < 1mm                               28.939     10.003   2.893
## logdose.particles.mL.master                      2.912      1.001   2.908
## logdose.mg.L.master:size_f100nm < 1µm            3.637      1.165   3.123
## logdose.mg.L.master:size_f1µm < 100µm            3.002      1.078   2.785
## logdose.mg.L.master:size_f100µm < 1mm            2.298      1.148   2.002
## size_f100nm < 1µm:logdose.particles.mL.master   -2.905      1.047  -2.773
## size_f1µm < 100µm:logdose.particles.mL.master   -2.767      1.004  -2.756
## size_f100µm < 1mm:logdose.particles.mL.master   -1.988      1.153  -1.725
##                                               Pr(>|z|)   
## (Intercept)                                    0.00227 **
## logdose.mg.L.master                            0.00672 **
## size_f100nm < 1µm                              0.00453 **
## size_f1µm < 100µm                              0.00390 **
## size_f100µm < 1mm                              0.00382 **
## logdose.particles.mL.master                    0.00364 **
## logdose.mg.L.master:size_f100nm < 1µm          0.00179 **
## logdose.mg.L.master:size_f1µm < 100µm          0.00535 **
## logdose.mg.L.master:size_f100µm < 1mm          0.04526 * 
## size_f100nm < 1µm:logdose.particles.mL.master  0.00555 **
## size_f1µm < 100µm:logdose.particles.mL.master  0.00585 **
## size_f100µm < 1mm:logdose.particles.mL.master  0.08459 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 945.92  on 853  degrees of freedom
## Residual deviance: 904.63  on 842  degrees of freedom
## AIC: 928.63
## 
## Number of Fisher Scoring iterations: 5

The above glm is visualized below.

acute_discrete <- ggpredict(m1_crust_acute_model_discrete,terms = c(#"logdose.mg.L.master", 
                                                                    "logdose.particles.mL.master",
                                                                    "size_f")) #ensure discrete factor is second
plot(acute_discrete)#, facet = TRUE, add.data = TRUE)

3.4.0.2 Alt Method using GG Plot

This approach achieves a similiar product as using the ggPredict() function, except it relies on ggplot(), which is more malleable and transparent. The general steps are to first create a new dataframe over 1000 values of size using expand.grid() then use predict() and plot() with geom_line() and colour=size.

require(viridis)
#filtered dataset
#m1_crust
#model
#summary(m1_crust_model_dose )
#generate distribution of data
mockData <- expand.grid(size.length.um.used.for.conversions = seq(0.034, 5000, 1),
            logdose.particles.mL.master = seq(-4.195, 12.650,0.1))

mockData$effect_10 <- predict(m1_crust_model_dose,
                     mockData,
                     type = "response")
#plot
mockData %>% 
  ggplot(aes(x = logdose.particles.mL.master, 
             y = effect_10, 
             color = size.length.um.used.for.conversions)) +
  geom_line() +
  scale_color_gradientn(colors = topo.colors(7)) +
  #scale_color_viridis_c(option = "A") +
  dark_theme_bw()

3.5 Survival package

3.5.1 Cumulative Hazard by Time

#include time
m1_crust <- aoc_setup %>% 
  filter(!effect_f == "NA") %>% 
  #filter(effect.metric == c("HONEC", "LOEC", "NOEC")) %>% 
  #filter(organism.group == "Crustacea") %>% 
  filter(!acute.chronic_f == "Unavailable") %>% 
  filter(bio.org == "organism") %>% 
  filter(!environment == "Terrestrial") %>% 
  filter(lvl1_f == "Fitness") %>% 
  filter(!size_f == "Not Reported") %>% 
  mutate(logdose.mg.L.master = log10(dose.mg.L.master)) %>% 
 # filter(acute.chronic_f == "Acute") %>% 
  dplyr::select(c(effect_10, logdose.mg.L.master, size_f, exposure.duration.d)) %>% 
  drop_na()

require(survival)
survival <- survival::coxph(survival::Surv(exposure.duration.d, effect_10) ~ logdose.mg.L.master + size_f, data = m1_crust)
#cumulative hazard
cumHaz <- ggpredict(survival,terms = c("logdose.mg.L.master", "size_f"), type = "cumhaz")
plot(cumHaz, facet = TRUE, add.data = TRUE, colors = "flat")

### Probability of Survival by Time

#survival probability over time
pr <- ggpredict(survival, c("logdose.mg.L.master", "size_f"), type = "surv")
plot(pr, colors = "social")